#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>

using namespace std;

class imagefile {
  int sizex;
  int sizey;
  int maxval;
  ifstream fin;
  unsigned char *r, *g, *b;
public:
  imagefile(const char *filename);
  unsigned char rpix(int x, int y);
  unsigned char gpix(int x, int y);
  unsigned char bpix(int x, int y);
  int sx();
  int sy();
  ~imagefile();
};


class imageiter {
  imagefile &F;
  int iterx, itery;
public:
  imageiter(imagefile &_F);
  int px();
  int py();
  void advance();
  bool eoi();
};

imagefile::imagefile(const char *filename) :
//  sizex(0), sizey(0), fin(filename,ios::binary|ios::in) {
  sizex(0), sizey(0), fin(filename) {
  if (!fin.is_open()) {
    cerr << "*** ERROR: Could not open input file " << filename << "\n";
    exit(0);
  }
  char line[100];
  fin.getline(line, 100);
  if ((line[0] != 'P') || (line[1] != '6')) {
    cerr << "*** Input is not a raw PPM file\n";
    fin.close();
    exit(0);
  }
  char c;
  char flag_c;
  while ((c = fin.get()) == '#') { //--- comment
    fin.getline(line, 100);
  }
  int maxval;
  fin.seekg(-1, ios::cur);
  fin >> sizex >> sizey >> maxval;
  fin.get(); // last NL
  cerr << "*** INFO: PPM file " << filename;
  cerr << " contains " << sizex << " by " << sizey << " pixels";
  cerr << " with range " << maxval << "\n";
  r = new unsigned char[sizex * sizey];
  g = new unsigned char[sizex * sizey];
  b = new unsigned char[sizex * sizey];
  for (int i=0; i<sizex*sizey; ++i) {
    fin.unsetf(ios::skipws); fin >> r[i] >> g[i] >> b[i];
  }
  fin.close();
}

unsigned char imagefile::rpix(int x, int y) {
  if ((x>=0) && (x < sizex) && (y >= 0) && (y < sizey))
    return r[y*sizex + x];
  else
    return 0;
}

unsigned char imagefile::gpix(int x, int y) {
  if ((x>=0) && (x < sizex) && (y >= 0) && (y < sizey))
    return g[y*sizex + x];
  else
    return 0;
}

unsigned char imagefile::bpix(int x, int y) {
  if ((x>=0) && (x < sizex) && (y >= 0) && (y < sizey))
    return b[y*sizex + x];
  else
    return 0;
}  
  
int imagefile::sx() {
    return sizex;
}

int imagefile::sy() {
  return sizey;
}

imagefile::~imagefile() {
  if (r) delete r;
  if (g) delete g;
  if (b) delete b;
}



imageiter::imageiter(imagefile &_F) : F(_F) {
  iterx = 0;
  itery = 0;
}

int imageiter::px() {
  return iterx;
}
 
int imageiter::py() {
  return itery;
}

void imageiter::advance() {
  ++iterx;
  if (!(iterx % 8)) {
    // crossed edge
    iterx -= 8;
    
    ++itery;
    if (!(itery % 8)) {
      // crossed edge - block complete
      itery -= 8;
      
      // advance block
      iterx += 8;
      if (iterx >= F.sx()) {
	// line wrap
	iterx = 0;
	
	itery += 8;
	if (itery >= F.sy()) {
	  // image complete
	  itery = 0;
	}
      }
    }
  }
}

bool imageiter::eoi() {
  return ((iterx == 0) && (itery == 0));
}


//-----------------------------------------------------------------

int main(int argc, char *argv[]) {
  
  if (argc != 3) {
    cerr << "Usage: " << argv[0] << " <PPM_file1> <PPM_file2>\n";
    exit(0);
  }

  imagefile F1(argv[1]);
  imagefile F2(argv[2]);

  if ((F1.sx() != F2.sx()) || (F1.sy() != F2.sy())) {
    cerr << "IMages have different size - cannot compare\n";
    exit(0);
  }

  imageiter I1(F1);
  imageiter I2(F2);

  double e = 0.;
  double re, ge, be;

  do {
    re = (F1.rpix(I1.px(), I1.py()) - F2.rpix(I2.px(), I2.py()));
    ge = (F1.gpix(I1.px(), I1.py()) - F2.gpix(I2.px(), I2.py()));
    be = (F1.bpix(I1.px(), I1.py()) - F2.bpix(I2.px(), I2.py()));
    I1.advance();
    I2.advance();

    e += re*re + ge*ge + be*be;
  } while (! I1.eoi());

  // normalize
  e = e / (255. * 255. * 255. * F1.sx() * F1.sy());

  double le;

  if (e < 1e-10) {
    cout << "SNR > 150 dB\n";
  } else {
    le = 10 * log10(e);
    cout << "SNR (dB) " << -le << "\n";
  }

  return 0;
}


